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ABSTRACT 

Analytic solutions for the evolution of cosmological linear density perturbations in the 
baryonic gas and coUisionless dark matter are derived. The solutions are expressed in a 
closed form in terms of elementary functions, for arbitrary baryonic mass fraction. They 
are obtained assuming f2 = 1 and a time independent comoving Jeans wavenumber, kj . 
By working with a time variable r = \a{t^/^), the evolution of the perturbations is 
described by linear differential equations with constant coefficients. The new equations 
are then solved by means of Laplace transformation assuming that the gas and dark 
matter trace the same density field before a sudden heating epoch. In a dark matter 
dominated Universe, the ratio of baryonic to dark matter density perturbation decays 
with time roughly like exp(— 5r/4) cx to the limiting value 1/[1 + (fc/fcj)^]. For 

wavenumbers k > kj / \/24, the decay is accompanied with oscillations of a period 



B7r/v/24(fc/fc,7 



1 in T. In comparison, as r increases in a baryonic matter dominated 



Universe, the ratio approaches 1 — (fc/fcj)^ for k < kj, and zero otherwise. 
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1 INTRODUCTION 

Several methods (e.g., Croft et al. 1998; Gnedin 1998; Nusser & Haehnelt 1999) have recently been proposed for extracting 
information on the mass density field from the Lyman-a forest. The underlying physical picture behind these method is 
that the absorbing Neutral hydrogen in the low density intergalactic medium (IGM) is tightly related to the mass density 
field on scales larger than the Jeans length. Below the Jeans length gas pressure segregates the baryons from the total mass 
fluctuations. On scales near the Jeans length, the evolution of the baryonic perturbation can affect estimates of the clustering 
amplitude from observations of the Lyman-a forest. 

Hydrodynamical simulations (Petitjean, Miicket & Kates 1995; Zhang, Anninos & Norman 1995; Hernquist et al. 1996; 
Miralda-Escude et al. 1996, Theuns et al. 1998) and semi-analytic models (e.g., Bi et. al. 1992, Gnedin & Hui 1996) of the 
IGM have been successful at explaining observations of the forest. Despite the success of the simulations, it is usually difficult 
to use them to study in detail the evolution of the gas below Jeans length (Theuns et al. 1998). The equations governing the 
evolution of baryons and dark matter in the nonlinear regime are extremely difficult to solve, even for special configurations 
like spherical collapse. Fortunately, since most of the IGM is of moderate density, linear analysis can be a suitable tool for 
understanding the evolution of the baryons (Gnedin & Hui 1998). Here we derive analytic solutions to the linear equations 
in a flat universe without a cosmological constant. Although the linear equations can readily be numerically integrated under 
a variety of conditions (Gnedin & Hui 1998), analytic treatment offers better understanding of the equations. Further, the 
paucity of analytic solutions makes their pursue worthwhile, even if tedious at times. The analytic solutions we derive here 
are subject to the condition that the baryonic and dark matter trace the same density and velocity fields before a sudden 
reionization epoch. After reionization the temperature of the IGM is assumed to be inversely proportional to the scale factor 
so that the comoving Jeans length is constant. 

The paper is organized as follows. In section 2 we cast the equations in the form of linear differential equations with 
constant coefficients. In Section 3 we present the solutions to these equations for several cases. We conclude in Section 4. 
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2 THE LINEAR EQUATIONS 

Let 5x{t, k) and 5i,{t, k) be , respectively, the Fourier modes of baryonic and dark matter density fluctuations. Let also fx and 
fb = 1 — fx he the the mean mass fractions of these two types of matter. We will restrict the analysis to perturbations in a 
flat universe without a cosmological constant. The linear equations governing the evolution of 5b and 5x are (e.g., Bi et. al. 
1992, Padmanabhan 1993, Gnedin & Hui 1998), 

where a[t) oc t^^^ is the scale factor, H{t) = 2/(3f) the Hubble function, and kj is the comoving Jeans wavenumber related 
to the speed of sound and the mean total (baryons plus dark) density, p = 3H^ /{8nG), by 

f3aH 



For fx close to unity the gravity of the baryons is negligible and Sx oc a. If also kj — const (i.e., oc 1/a) and we impose, 
at some initial time ti, the condition Sb{ti) — 5x{ti)/{l + {k/kj)"^), then 5b{t) — 5xit)/{l + (k/kj)^) at any later time, t. In 
the next section we will show how to solve these equations with < lefx < assuming that Sb{t) — Sx{t) for t < ti, which is 
appropriate for a sudden reionization of the IGM at ti. 

In order to solve these equations we work with a new time variable r = ln(a) instead of t (Nusser & Colberg 1998). In 
terms of r the equations become, 

d^Sx ^IdSx 3 

d'^Sb 1 dfJb 3 , , 3 2 c 

d7^+2d7 ^ ^if^Sx + fb5b)^-nSb, (3) 
where have defined k = k/kj. 

When hi is constant, the differential equations ^ are linear with constant coefficients and they can be solved by means of 
Laplace transformation. Since Laplace transforms are seldom used in cosmological studies, it seems prudent to briefiy review 
their basic properties which are relevant to us. We refer the reader to Arfkin (1985) and references therein for mathematical 
details. The Laplace transform, /(s), of a function F(t), where t > 0, is defined as 

/•oo 

f{s)=C{F{t)}^ / exp(-st)F(t)dt (4) 
Jo 

We will need the Laplace transforms of first and second derivatives of a function. Using (^ these transforms can be related 
to the /(s) by 

aF'{t)} = sf{s)-F{Q) 

C{F"(t)} = sV(s) - sF(0) - F'(0) , (5) 

where the prime and double prime denote first and second order derivatives, respectively. The Bromwich integral expresses 
Fit) in terms of f[s) as 

nt)=2;-/ exp(st)/(s)ds (6) 



7 — lOO 



where i — \/— T and 7 is a real number chosen so that all poles of /(s) lie, in the complex plane, to the left of the vertical line 
defining the integration path. Therefore, by the residue theorem we have 

F{t) = ^ [residues of exp (st) f{s)] (7) 

As an example consider f{s) — l/(s — si)(s — S2) which has two simple poles at s = si and S2. The residues of exp{st)f{s) at 
these poles are exp(sif) /(si — S2) and — exp(s2t)/(si — S2) so that, by (^), F{t) = [exp(sit) — exp(s2f)]/(si — S2). If si = S2, the 
function has a pole of order two at si. The residue in this case is d[(s — si)^ exp(st)/(s))]/ds evaluated at s = si. Therefore, 
F{t) = texp(sit). 
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3 THE SOLUTIONS 

Denote by and A^ the Laplace transforms of 5x and St, respectively. By taking the Laplace transform of we can 
obtain relations between A^ and A;,. The initial conditions are contained in the Laplace transforms of the first and second 
derivatives of the densities. So first we have to specify in mathematical terms our choice for the initial conditions. For simplicity 
of notation we fix the initial conditions at r = 1 assuming that before that time the temperature of the baryonic fiuid is zero, 
i.e., K — 0. The initial conditions are fixed by the values of 5^ and 5b and their first derivatives at r = 1. Before r = 1, we 
have 5x = Sb = exp(T), ignoring the decaying mode and setting arbitrarily 5x{t — 1) = 1. The first derivatives of Sx and 5b are 
therefore equal to unity at r = 1. This fixes the initial conditions necessary for solving (^). Although we will present solutions 
satisfying only these initial conditions, we will, for completeness, write the Laplace transformation of for 5b = aS^ and 
d(5(,/dr = ad5x/dT, at r = 1. Then the Laplace transformation of yields 

(s^' + is^A, = ^ (/,A, + /bAf,) + s + I 

(s^ + is^Ai = ^ifxAx+ fbAb)-^K^Ab + as+^a , (8) 

where have have used to computed the transforms of the first and second derivatives of 5x and 5b. For = 1 the first of 
these equations yields 



1 



(9) 



S — 1 

which is the Laplace transform of exp(r). If we take a = {1 + k'^)"^ and substitute ^ in the second equation of we get 
1 1 A, 



S - 11 + k2 1 + k2 ' 

,2\ 



(10) 



which leads to the well known solution 5b — Sx / {1 + n ). 

Subsequently we will present solutions only for a = 1. In this case, equations (S) yield 



(11) 



and 

Before solving these equations for any value of fx in the range 0-1, it is instructive to examine the solutions for the 
special values fx = l and 0. 



3.1 Case I: fx = 1 

In this case A^, = (s — 1)~^ and equation ( p^ ) can be written in the form 



(13) 



(s-l)(s-s_)(s-s^) 
where are the roots of /2 + + s/2. They are given by 

s±=-i(l±x) ; x' = 1-24^2 (14) 

We will deal with the case = at the end of this subsection. For 7^ 0, all three poles of Aj, are simple and so, by (^, 
its inverse transform is 



exp (r) 1 
06 = — - + 



1 + S_ — 1 + 



(s_ — l) exp {s^r'j — (s^ — l) exp (•s_''") 



(15) 
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For > the roots are real and the result is 



exp (r) 1 
Ob = — : — — + 



1 + k2 



2x1 + ^^ 



(x-5)exp(-|r) +(5 + x)exp(+Jr) exp(-jj 



(16) 



The first term is the solution given in the previous section for a — 1/(1 + k^). The maximum value attains is unity, so the 
second and third terms are always decaying. If x^ < 0, then ( [l5| ) gives the solution 



exp(r) 1 k"^ 
Oh = — \ — + 



5sin(jr)+xcos(jr) exp (-l) 



(17) 



where x is the imaginary part of x- The solution shows an oscillatory behavior with a period of IGtt/x- The envelope of these 
oscillations decays like exp(— r/4) oc t~^^^ . 



We deal now with the case x^ = 0, which occurs for — 1/24. Here special care is needed because s_ 



-1/4. 



However the contribution of the pole at s = — 1 to the Bromwich integral remains unchanged and the contribution of the 
second order pole at s = —1/4 is simply the first derivative of exp(sr)(s + 1/4)^A6 at s = —1/4. The result is 



24 

= ^ exp (r) 



1 

20 



exp 



(18) 



The first term on the left is the familiar Sx/i'i + i<^) evaluated at — 1/24. The expression can also be derived by taking 
the limit x^ ^ in either (lli) or (111). 



In the limit r oo the ratio 5i,/5x is 1/(1 + k'^). 



3.2 Case II: ^ = 

This is equivalent to ignoring the gravity of the dark matter. Of course here only the behavior of the perturbation in the 
baryons is relevant since the dark matter plays no role. However, for the sake of completeness and comparison with other 
situations we will solve for the dark matter fluctuations as well. 

We first find the solution for Sx- If fx ~ 0, we can express ( pT| ) in terms of , the roots of Zk^ /2 + + s/2 — 3/2, as 



{s-S_) {S-S^) 



(19) 



where 



5± =-4(l±x) 



X^ = 25 - 24k^ 



(20) 



If K 7^ 1, then the function A;, has three simple poles at s = 0, s_, and s^. So for k 7^ 1 and x^ > 0, Bromwich integral 
yields 



2exp(-0 



^ 1-^2 

The expression for Si, is 



-3 



+ 



2x 1 - 



(X - 5) exp + (x + 5) exp (+f exp (-0 



(21) 



2X 



(x-5)exp(-Jr) + (x + 5) exp (|r) exp(-l) 



(22) 



When K = 1 the solution can be found either by taking the limit k ^ 1 in the previous two expressions or by direct 
evaluation of the Bromwich integral with a two poles of order two at s = and —1/2. The result is 



Sx = -27 + 28exp ( -- ) +3r 



3 + 2 exp - 



and 



Sh — i ~ 2 exp ( — 



(23) 
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This implies that Sx grows hnearly with r at late times while 5b reaches an asymptotic value of 3. 

The oscillatory behavior of St and Sx appears when < 0, i.e., for > 25/24. The expressions in this case can be 
obtained by replacing x with ix in (^l|) and (|22|). For » 25/24, the solution coincides with that given in Padmanabhan 
(1993) 

In the limit t oo the ratio St/Sx is 1 — for k < 1 and, zero otherwise. 



3.3 Case III: < fx < 1 

Again we first derive Sx- The denominator and numerator in ( [ill ) do not have any common roots. Then the poles of Ax are 
the roots of the denominator. We find these roots as follows. Denote y{s) = + s/2 — 3/2 and equate denominator to zero 
to obtain 

y^ + ^y(l + ^^)+^ft^=0 (24) 

where we have used 1 — fx = fb- This equation is satisfied for the following values of y, 

y,,^ = -j{l + ^){l±E) ; H^-1-4/. ^\ , (25) 

4 ^ ' (1 + K'') 

where the subscripts p and m correspond to the plus and minus sign, respectively. So the roots of the denominator in ( pi] ) 
are the values of s which make y(s) = yp,m- Let Sp^± and ^ be the roots of y{s) ~ yp — and y{s) — ym ~ 0, respectively. 
They are given by 

s,„,± =-i(l±Xm) ; Xm' = 25-12(l + A^2)(l + H) , 
Sp,± =-j(l±Xp) ; Xp' = 25-12(l + ,c^)(l-H) . (26) 

Excluding the values fb = and 1, which we have considered in the previous subsections, we have < < 1 - /i,. This 
ensures that all four roots, ^ and ^ , are distinct. Also, 1 < Xp < 25 for any k, so the roots ^ are real. On the other 
hand, < Xm < 1 when < 25/(600 — 576fb), and negative otherwise. So ^ can be complex. 

For Xm > 0, the Bromwich integral yields 



^ _ 24k^+xI-1 
2xp (Xp - Xm) 



(Xp-5)exp(-^r) +(xp + 5)exp(+^r) exp (-0 + (xp ^ Xm) , (27) 



with the second term on the r.h.s is obtained interchanging Xp a^nd Xm in the first term. This result can be extended to 
Xm < by writing Xm = iXm where Xm is real. 

Using (^2|) we similarly obtain 5b 

(Xp-5)exp(-^r) +(xp + 5)exp(+^T) exp (-l) + (xp ^ Xm) . (28) 



x^ 



2Xp (Xp - Xm) _ 

To visualize these solutions, we show in Fig.l the density evolution for various values of fx and k. The dark matter curve 
for fx = 0.9 is very close to exp(r). For k = 5 the baryonic perturbations for both values of fx show oscillations with similar 
period and amplitude. This is simply because for high k, the evolution is mainly dictated by pressure forces. 

In Fig. 2 we the solid lines represent the ratio 5b{n) / 5x{n) as a function of k at different t designated in the plot by the 
redshift, z. The initial conditions are satisfied at z = 6 and fx ~ 0.9 was taken. Also plotted, as the dotted line in each panel, 
the function 1/(1 + k^) which represents the limiting solution as r ^ oo. The analytic curves show more oscillations as they 
get closer to the limiting solution, when the redshift is decreased. 

In many applications (e.g., Bi et al 1992, Bi & Davidsen 1997, Nusser & Haehnelt 2000) the limiting ratio 1/(1 + k^) is 
often used to filter the mass power spectrum in order to generate density fiuctuations in the gas. As pointed out by Gnedin 
& Hui (1998) this may lead to a significant bias in statistics of the gas density. As an illustration of this bias we compute the 
rms values of the gas density by filtering a scale free mass power spectrum of slope n — —2.5 with the ratio {5b/5x)'^ given 
from the analytic, and the limiting solutions, respectively. Fig. 3 shows the ratio of the former to the latter rms value. As in 
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the previous figure, fx = 0.9 and initial conditions satisfied at 2 = 6. Using the filter 1/(1 + k^) can seriously underestimate 
the amplitude of gas density fluctuations. Only when we approach z = the ratio gets close to to unity. 



4 SUMMARY 

We have found analytic solutions to the linear equations governing the evolution of baryonic and dark matter under four 
assumptions. First, the Universe is flat without cosmological constant. Second, sudden reionization of the IGM. Third, the 
temperature of the low density IGM drops like 1/a so that the comoving Jeans length is time independent. Fourth, before 
reionization the IGM is cold and the baryonic and dark matter trace the same density and velocity fields. 

Of these assumptions, only the fourth has a physical basis, at least before any heating has occurred and when the IGM 
temperature is low. This is also the only assumption which if changed, the equations can still be readily solved by Laplace 
transformation. Unfortunately, relaxing any of the other assumptions complicates the analytic treatment of equations by 
means of Laplace transformation. For example, suppose that the Jeans wavenumber changes with time according to . Then 
the Laplace transformation of the term involving kj will yield A;, at (s + (3) while other terms involve Ai,(s). 

Yet the solutions can be useful for semi-analytic modeling of the IGM. They offer a convenient improvement over the 
commonly used filter 1/[1 + (fe/fej)^] for generating gas fluctuations associated with a given mass density field. 

The analytic solutions presented here were verified by a comparison with the solutions obtained by numerical integration 
of equations Aw. All numerical and analytic solutions agreed up to the numerical accuracy. 
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Figure 1. Curves of Sx{t) (solid lines) and Si,{t) (dotted) for various values of fx and k = k/kj. 
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Figure 2. Ratios of baryonic to dark matter density as a funetion of k at different times, indicates by the redshift, z, in each paneL 
The sohd lines arc the analytic solutions obtained with fx = 0.9 and where the initial conditions are satisfied at z = 6. The dotted line 
in each panel shows 1/(1 + k^), the ratio corresponding to the limiting solution. 
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Figure 3. The ratio of tlie baryonic density rms values computed by filtering the mass power spectrum with (<5i,(T, K)/6x{t, n))^ in 
analytic solution to that obtained by filtering with 1/(1 + /t^)^. A scale free mass power spectrum with slope n = —2.5 is assumed. 



